Lab 5

Author

Giuliet Kibler

Learning goals

  • Use the merge() function to join two datasets.
  • Deal with missings and impute data.
  • Identify relevant observations using quantile().
  • Practice your GitHub skills.

Lab description

For this lab we will be, again, dealing with the meteorological dataset downloaded from the NOAA, the met. In this case, we will use data.table to answer some questions regarding the met dataset, while at the same time practice your Git+GitHub skills for this project.

This markdown document should be rendered using github_document document.

Part 1: Setup the Git project and the GitHub repository

  1. Go to your documents (or wherever you are planning to store the data) in your computer, and create a folder for this project, for example, “PM566-labs”

  2. In that folder, save this template as “README.Rmd”. This will be the markdown file where all the magic will happen.

  3. Go to your GitHub account and create a new repository, hopefully of the same name that this folder has, i.e., “PM566-labs”.

  4. Initialize the Git project, add the “README.Rmd” file, and make your first commit.

  5. Add the repo you just created on GitHub.com to the list of remotes, and push your commit to origin while setting the upstream.

Most of the steps can be done using command line:

# Step 1
cd ~/Documents
mkdir PM566-labs
cd PM566-labs

# Step 2
wget https://raw.githubusercontent.com/USCbiostats/PM566/master/website/content/assignment/05-lab.Rmd 
mv 05-lab.Rmd README.md

# Step 3
# Happens on github

# Step 4
git init
git add README.Rmd
git commit -m "First commit"

# Step 5
git remote add origin git@github.com:[username]/PM566-labs
git push -u origin master

You can also complete the steps in R (replace with your paths/username when needed)

# Step 1
setwd("~/Documents")
dir.create("PM566-labs")
setwd("PM566-labs")

# Step 2
download.file(
  "https://raw.githubusercontent.com/USCbiostats/PM566/master/website/content/assignment/05-lab.Rmd",
  destfile = "README.Rmd"
  )

# Step 3: Happens on Github

# Step 4
system("git init && git add README.Rmd")
system('git commit -m "First commit"')

# Step 5
system("git remote add origin git@github.com:[username]/PM566-labs")
system("git push -u origin master")

Once you are done setting up the project, you can now start working with the MET data.

Setup in R

  1. Load the data.table (and the dtplyr and dplyr packages if you plan to work with those).

    library(data.table)
    library(dtplyr)
    library(dplyr)
    
    Attaching package: 'dplyr'
    The following objects are masked from 'package:data.table':
    
        between, first, last
    The following objects are masked from 'package:stats':
    
        filter, lag
    The following objects are masked from 'package:base':
    
        intersect, setdiff, setequal, union
  2. Load the met data from https://raw.githubusercontent.com/USCbiostats/data-science-data/master/02_met/met_all.gz, and also the station data. For the later, you can use the code we used during lecture to pre-process the stations data:

    # Lead met data
    if (!file.exists("met_all.gz"))
      download.file(
        url = "https://raw.githubusercontent.com/USCbiostats/data-science-data/master/02_met/met_all.gz",
        destfile = "met_all.gz",
        method   = "libcurl",
        timeout  = 60
        )
    met <- data.table::fread("met_all.gz")
# Download the data
stations <- fread("ftp://ftp.ncdc.noaa.gov/pub/data/noaa/isd-history.csv")
stations[, USAF := as.integer(USAF)]
Warning in eval(jsub, SDenv, parent.frame()): NAs introduced by coercion
# Dealing with NAs and 999999
stations[, USAF   := fifelse(USAF == 999999, NA_integer_, USAF)]
stations[, CTRY   := fifelse(CTRY == "", NA_character_, CTRY)]
stations[, STATE  := fifelse(STATE == "", NA_character_, STATE)]

# Selecting the three relevant columns, and keeping unique records
stations <- unique(stations[, list(USAF, CTRY, STATE)])

# Dropping NAs
stations <- stations[!is.na(USAF)]

# Removing duplicates
stations[, n := 1:.N, by = .(USAF)]
stations <- stations[n == 1,][, n := NULL]
  1. Merge the data as we did during the lecture.

    # Convert the USAF column to integer
    stations$USAF <- as.integer(stations$USAF)
    # Check for and delete duplicates in stations
    stations[, n := 1:.N, by = .(USAF)]
    stations <- stations[n == 1,][, n := NULL]
    
    # Merge the data
    dat <- merge(
      # Data
      x     = met,      
      y     = stations, 
      # List of variables to match
      by.x  = "USAFID",
      by.y  = "USAF", 
      # Which obs to keep?
      all.x = TRUE,      
      all.y = FALSE
      )

Question 1: Representative station for the US

What is the median station in terms of temperature, wind speed, and atmospheric pressure? Look for the three weather stations that best represent continental US using the quantile() function. Do these three coincide?

# Calculate medidians by station
medians_by_USAFID <- aggregate(
  cbind(temp, wind.sp, atm.press) ~ USAFID, 
  data = dat, 
  FUN = median, 
  na.rm = TRUE
)
# Calculate median of the country 
median_temp <- quantile(dat$temp, 0.5, na.rm = TRUE)
median_wind <- quantile(dat$wind.sp, 0.5, na.rm = TRUE)
median_press <- quantile(dat$atm.press, 0.5, na.rm = TRUE)

# Find the median stations
closest_temp <- medians_by_USAFID[which.min(abs(medians_by_USAFID$temp - median_temp)), ]
closest_wind <- medians_by_USAFID[which.min(abs(medians_by_USAFID$wind.sp - median_wind)), ]
closest_press <- medians_by_USAFID[which.min(abs(medians_by_USAFID$atm.press - median_press)), ]

closest_temp
    USAFID  temp wind.sp atm.press
229 722860 23.55     1.5    1012.5
closest_wind
  USAFID temp wind.sp atm.press
5 720333 22.8     2.1    1012.4
closest_press
    USAFID temp wind.sp atm.press
131 722420 30.6     4.6    1014.1

The median station for temp is USAFID 722860 in California, for wind speed is USAFID 720333 in California, and for atmospheric pressure is USAFID 722420 in Texas. None of these medians coincide.

Knit the document, commit your changes, and push it to GitHub. Don’t forget to add README.md to the tree, the first time you render it.

Question 2: Representative station per state

Just like the previous question, you are asked to identify what is the most representative, the median, station per state. This time, instead of looking at one variable at a time, look at the euclidean distance. If multiple stations show in the median, select the one located at the lowest latitude.

# Find the medians of each state
medians_by_state <- dat |>
  group_by(STATE) |>
  summarize(
    median_temp = median(temp, na.rm = TRUE),
    median_wind = median(wind.sp, na.rm = TRUE),
    median_press = median(atm.press, na.rm = TRUE)
  )
# Find the median stations
closest_usafid_by_state <- dat |>
  group_by(STATE)|>
  reframe(
    closest_usafid_temp = USAFID[which.min(abs(temp - medians_by_state$median_temp[match(STATE, medians_by_state$STATE)]))],
    closest_usafid_wind = USAFID[which.min(abs(wind.sp - medians_by_state$median_wind[match(STATE, medians_by_state$STATE)]))],
    closest_usafid_press = USAFID[which.min(abs(atm.press - medians_by_state$median_press[match(STATE, medians_by_state$STATE)]))]
  )
closest_usafid_by_state
# A tibble: 46 × 4
   STATE closest_usafid_temp closest_usafid_wind closest_usafid_press
   <chr>               <int>               <int>                <int>
 1 AL                 720265              720265               720361
 2 AR                 720175              720172               720175
 3 AZ                 720339              720339               722720
 4 CA                 690150              690150               690150
 5 CO                 720528              720385               722817
 6 CT                 725027              720545               725027
 7 DE                 724088              724088               724093
 8 FL                 720373              720373               720383
 9 GA                 720257              720257               722070
10 IA                 720293              720293               725420
# ℹ 36 more rows

Knit the doc and save it on GitHub.

Question 3: In the middle?

For each state, identify what is the station that is closest to the mid-point of the state. Combining these with the stations you identified in the previous question, use leaflet() to visualize all ~100 points in the same figure, applying different colors for those identified in this question.

# Calculate the midpoints for each state
state_midpoints <- dat |>
  group_by(STATE) |>
  summarize(
    mid_lat = mean(lat, na.rm = TRUE),  # Replace with your latitude column name
    mid_lon = mean(lon, na.rm = TRUE), # Replace with your longitude column name
    .groups = 'drop'
  )

# Find the closest station to each state midpoint
closest_to_midpoint <- state_midpoints |>
  rowwise() |>
  mutate(
    USAFID = dat$USAFID[which.min(sqrt((dat$lat - mid_lat)^2 + (dat$lon - mid_lon)^2))]
  ) |>
  ungroup()

# Combine the closest stations for graphing
combined_stations <- closest_to_midpoint |>
  left_join(dat |> select(USAFID, lat, lon), by = "USAFID")

# Graph using leaflet()
library(leaflet)

# Create a leaflet map
map <- leaflet() |>
  addTiles()  # Add OpenStreetMap tiles

# Add points for closest stations to midpoints
map <- map |>
  addCircleMarkers(data = combined_stations,
                   ~lon, ~lat,
                   color = "red",  # Color for midpoints
                   radius = 5,
                   popup = ~paste("Station ID:", USAFID))
map <- map |>
  addCircleMarkers(data = state_midpoints,
                   ~mid_lon, ~mid_lat,
                   color = "blue",  # Color for state midpoints
                   radius = 7,
                   label = ~paste("State Midpoint:", STATE),
                   group = "State Midpoints")
map

Knit the doc and save it on GitHub.

Question 4: Means of means

Using the quantile() function, generate a summary table that shows the number of states included, average temperature, wind-speed, and atmospheric pressure by the variable “average temperature level,” which you’ll need to create.

Start by computing the states’ average temperature. Use that measurement to classify them according to the following criteria:

  • low: temp < 20
  • Mid: temp >= 20 and temp < 25
  • High: temp >= 25
# Averages by state
means_by_state <- dat |>
  group_by(STATE) |>
  summarize(
    mean_temp = mean(temp, na.rm = TRUE),
    mean_wind = mean(wind.sp, na.rm = TRUE),
    mean_press = median(atm.press, na.rm = TRUE)
  )

# Categorize by average temperature
means_by_state <- means_by_state |>
  mutate(avg_temp_level = case_when(
    mean_temp < 20 ~ "Low",
    mean_temp >= 20 & mean_temp < 25 ~ "Mid",
    mean_temp >= 25 ~ "High",
    TRUE ~ NA_character_  # This handles any NA values if needed
  ))

# Create a summary table by number of states in each temperature category
summary_table <- means_by_state |>
  group_by(avg_temp_level) |>
  summarize(
    number_of_states = n(),
    avg_temp = mean(mean_temp, na.rm = TRUE),
    avg_wind = mean(mean_wind, na.rm = TRUE),
    avg_press = mean(mean_press, na.rm = TRUE),
    .groups = 'drop'
  )
summary_table
# A tibble: 3 × 5
  avg_temp_level number_of_states avg_temp avg_wind avg_press
  <chr>                     <int>    <dbl>    <dbl>     <dbl>
1 High                         12     27.0     2.43     1014.
2 Low                          11     18.7     2.55     1014.
3 Mid                          25     22.6     2.39     1015.

Once you are done with that, you can compute the following:

  • Number of entries (records),
  • Number of NA entries,
  • Number of stations,
  • Number of states included, and
  • Mean temperature, wind-speed, and atmospheric pressure.

All by the levels described before.

# Categorize data by temperature 
dat <- dat |>
  mutate(avg_temp_level = case_when(
    temp < 20 ~ "Low",
    temp >= 20 & temp < 25 ~ "Mid",
    temp >= 25 ~ "High",
    TRUE ~ NA_character_  # This handles any NA values if needed
  ))

# Summarize by avg_temp_level
summary_stats <- dat |>
  group_by(avg_temp_level) |>
  summarize(
    number_of_entries = n(),
    number_of_na_entries = sum(is.na(dat)),
    number_of_stations = n_distinct(STATE), 
    mean_temperature = mean(temp, na.rm = TRUE),
    mean_wind_speed = mean(wind.sp, na.rm = TRUE),
    mean_pressure = mean(atm.press, na.rm = TRUE),
    .groups = 'drop'
  )
summary_stats
# A tibble: 4 × 7
  avg_temp_level number_of_entries number_of_na_entries number_of_stations
  <chr>                      <int>                <int>              <int>
1 High                      980808              2986380                 48
2 Low                       590757              2986380                 48
3 Mid                       745689              2986380                 48
4 <NA>                       60089              2986380                 48
# ℹ 3 more variables: mean_temperature <dbl>, mean_wind_speed <dbl>,
#   mean_pressure <dbl>

Knit the document, commit your changes, and push them to GitHub. Once you’re done, you can email me the link to your repository or, preferably, link to the Lab 5 Issue in your commit message, as follows:

git commit -a -m "Finalizing lab 5 https://github.com/USCbiostats/PM566/issues/68"

This will let me know which version of your repository to look at for grading purposes.